import numpy.fft as nf
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import PercentFormatter
import os

plt.rcParams['backend'] = 'PDF'
plt.rcParams['font.sans-serif'] = 'Arial Unicode MS'
plt.rcParams['axes.unicode_minus'] = False

path = '10Hz/'
files = os.listdir(path)

result = {}

for file in files:
    if file.startswith(".") or file.endswith("1") or file.endswith(".pdf") or file.endswith(".png") :
        continue
    data = np.loadtxt(f"{path}{file}", skiprows=9)  # 跳过前9行数据
    t = data[:, 0]  # 第一列数据为时间序列
    f_t = data[:, 1]  # 第二列数据为函数序列

    '''
    对自己定义的函数进行分析
    '''
    # B = 30.0    # 频谱分析的最大频率
    # f_s = 2 * B # 采样频率
    # delta_f = 0.01  # 频谱分辨率
    # N = int(f_s / delta_f)  # 采样点数
    # T = N / f_s # 采样时间

    # def signal_samples(t):
    #     # return (2 * np.sin(2*np.pi*t) + 3 * np.sin(22*2*np.pi*t) + 2 * np.random.randn(*np.shape(t)))
    #     return 2 * np.sin(4*2*np.pi*t) + 3 * np.sin(3*2*np.pi*t)

    # t = np.linspace(0, T, N)
    # f_t = signal_samples(t)

    fig, axes = plt.subplots(1, 2)  # 画两个子图，分别为时域图和频域图
    axes[0].plot(t, f_t)
    axes[0].set_xlabel("time (s)")
    axes[0].set_ylabel("signal")

    F = nf.fft(f_t)
    N = f_t.size  # 采样点数
    delta_t = t[1] - t[0]  # 采样间隔
    f = nf.fftfreq(f_t.size, delta_t)  # 频域的横坐标序列
    mask = np.where(f > 0)  # 只分析频率大于零的部分（小于零的部分与之完全对称）
    A = 2 * abs(F[mask]) / N  # 频谱图的纵坐标变换为实际振幅

    meas_freq = float(file.split('Hz')[0])  # 从文件名中读取激振器的频率值
    sign_freq = (f[f > 0])[np.argmax(abs(F[mask]))]  # 找到频谱图山振幅最大时对应的频率值
    amplitude = np.max(A)  # 最大振幅值
    print(f"分量最大的频率分量为：{sign_freq:.2f} Hz, 对应的幅值为：{amplitude:.2f} V")
    devi_freq = sign_freq - meas_freq  # 绝对误差
    real_devi = devi_freq / meas_freq  # 相对误差
    print(f"偏差为：{devi_freq:.2f} Hz（{real_devi:.2%}）")
    text = f"f = {sign_freq:.2f} Hz, A = {amplitude:.2f} V, error = {devi_freq:.2f} Hz, {real_devi:.2%}"
    plt.suptitle(text)  # 将误差作为图的标题，显示在图的上方

    axes[1].plot(f[mask], A)
    axes[1].set_xlabel("frequency (Hz)")
    axes[1].set_ylabel("$A$")
    plt.tight_layout(rect=[0, 0, 1, 0.92])  # 设置子图布局，避免重叠
    plt.savefig(f"{path}{meas_freq}Hz.pdf")  # 将分析结果保存在以频率信息命名的png图片中

    result.update({meas_freq: devi_freq})

li = list(zip(*sorted(result.items())))     # zip(*)相当于zip()的逆过程
li = np.array(li)

fig, axes1 = plt.subplots()
plt.title('测量结果误差分析')

color = 'tab:orange'
axes1.plot(li[0], li[1], color=color)
axes1.set_xlabel("频率 $f$/Hz")
axes1.set_ylabel("绝对误差 $\Delta$/Hz", color=color)
axes1.tick_params(axis='y', labelcolor=color)

axes2 = axes1.twinx()
color = 'tab:cyan'
axes2.plot(li[0], li[1] / li[0], color=color)
axes2.set_ylabel("相对误差 $\delta$", color=color)
axes2.yaxis.set_major_formatter(PercentFormatter(1))
axes2.tick_params(axis='y', labelcolor=color)

plt.tight_layout()
plt.savefig(f"{path}error_analyze.pdf")

plt.show()
